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Abstract What are the face-probabilities of a cuboidal die, i.e. a die with 
different side-lengths? This paper introduces a model for these probabilities 
based on a Gibbs distribution. Experimental data produced in this work and 
drawn from the literature support the Gibbs model. The experiments also 
reveal that the physical conditions, such as the quality of the surface onto 
which the dice are dropped, can affect the face-probabilities. In the Gibbs 
model, those variations are condensed in a single parameter, adjustable to the 
physical conditions. 
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1 Introduction 

Symmetry considerations make it seem obvious that an ideal cubical die lands 
on all six faces with an identical probability of g. However, what happens 
when non-fair dice are tossed? In particular, what are the face-probabilities 
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of a homogeneous cuboid, i.e. a six-sided polyhedron with parallel faces but 
different side-lengths? This is a surprisingly challenging problem. 

This paper offers a robust answer to the question above. It begins with a 
brief historical account (Section [2]) and presents a control experiment with a 
single cuboid (Section [3| that is used as a benchmark for the theoretical mod- 
elling. Section [4] then introduces a new model based on a Gibbs distribution, 
which is found to be consistent with the control experiment. This model natu- 
rally contains a free parameter, which characterizes the physical conditions of 
the experiment. The new model is then compared against experimental data 
with differently sized cuboids drawn from the literature (Section [5]) and ex- 
tended to non-cuboidal dice via the example of U-shaped dice (Section |6| . 
Section [7] summaries the findings of this paper. 



2 Brief history 

Already Isaac Newton mentioned the problem. In a private writing dated be- 
tween 1664 and 1666, and published in 1967 (Newton, 1967, p. 60-61), he 
wrote on the face-probabilities of a tossed cuboid: "if a die bee not a Regular 
body but a Parallelepipcdon or otherwise unequally sided, it may bee found 
how much one cast is more easily gotten then another." It remains unclear 
whether Newton really tried to solve this problem. 

In 1692, the problem appeared again in a paper by John Arbuthnot: "In a 
Parallelopipedon, whose Sides are to another in the Ratio of a,b,c: to find at 
how many Throws any one may undertake that any given Plane, viz. ab, may 
arise" (quoted from Hyksova et al., 2012). Arbuthnot wrote that he left "the 
solution to those who think it merits their pains." 

Fifty years later, Thomas Simpson (1740) used a simple geometrical idea 
to model the face-probabilities of a tossed cuboid. He assumed the probabil- 
ity of each face to be proportional to the surface area of the corresponding 
spherical quadrilateral, i.e. to the solid angle spanned by the face when seen 
from the centre of the cuboid. However, subsequent experimental investigations 
(e.g. Singmaster, 1981) clearly rejected Simpson's model. Budden (1980) and 
Heilbronner (1985) also experimented with series of cuboids. Although Budden 
and Heilbronner did not find a formula for the face-probabilities, their data 
again disqualifies the Simpson model, and so do modern computer simulations 
of tossed cuboids (Obreschkow, 2006). Regardless of the clear insufficiency of 
the Simpson model, a recent paper by Hyksova et al. (2012) still refers to this 
model without criticism. Because of this discrepancy, this paper will first reem- 
phasize the insufficiency of the Simpson model (Section [3]), before introducing 
a much more accurate model (Section [3]). 



3 Control experiment 

The control experiment is performed with a wooden (13 x 20 x 23 mm 3 )- 
cuboid. When tossing this cuboid, it became clear that the face-probabilities 
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Face i 


1 


2 


3 


4 


5 


6 


Surface area [mm 2 ] 


299 


260 


460 


460 


260 


299 


Half-height hi [mm] 


10 


11.5 


6.5 


6.5 


11.5 


10 


fi experiment I (N = 2, 700) [%] 


10.3 


7.7 


30.9 


32.7 


7.6 


10.9 


fi experiment II (N = 1,000) [%] 


5.5 


1.5 


43.5 


42.5 


2.6 


4.1 


Pi Simpson model [%] 


13.5 


10.5 


26.0 


26.0 


10.5 


13.5 


Pi Gibbs model (/3 = 4.90) [%] 


11.2 


7.2 


31.6 


31.6 


7.2 


11.2 


p t Gibbs model (fi = 10.2) [%] 


5.0 


2.0 


43.0 


43.0 


2.0 


5.0 



Table 1 Control experiment with a homogeneous (13 X 20 X 23 mm 3 )-cuboid. Faces 1 and 
6: 13 X 23 mm; faces 2 and 5: 13 X 20 mm, faces 3 and 4: 20 X 23 mm. See Section [3] for 
details. 



significantly depend on the physical conditions, such as the tossing technique, 
the height of free fall, the shape of the cuboid's edges, and the elasticity of 
the surface on which the cuboid lands. For example, a rough or elastic surface 
generally increases the face-probabilities of the two largest faces. The same 
qualitative change is observed when the cuboid is tossed from an arm-length 
above the table rather than using a dice cup. 

To account for the importance of the physical conditions, two experimental 
runs were performed. In experiment I, the cuboid was tossed N = 2, 700 times 
on a wooden table using a leather dice cup. In experiment II, the cuboid 
was dropped N = 1,000 times onto a polished steel surface from an initial 
height of 1 m. Table [I] lists the observed frequencies fi — rii/N, where i%i is 
the number of times that face i (i — 1,...,6) showed up. As expected, the 
measured frequencies differ significantly between experiment I and II, thus 
demonstrating that the shape of the cuboid alone does not determine the face- 
probabilities. Several physical reasons might be responsible for the different 
outcome probabilities in the two experiments. For example, the dice cup might 
have a stabilizing function when the cuboid lands on one of its small faces. In 
turn, dropping a cuboid from a high level (1 m) implies that the cuboid bounces 
off the floor many times before it comes to rest. Multiple bounces tend to result 
in a fast rolling around the longest axis of the cuboid, which implies that the 
cuboid is very unlikely to land on one of the two smallest faces. Qualitatively 
this suggests that the face-probabilities of the largest faces increase with the 
initial energy of the tossing process. 

The differences in the observed frequencies of opposite faces (e.g. face 3 
and 4) give a rough estimate of the deviation between measured frequencies 
and underlying probabilities. Those deviations would disappear if N — > oo. 

For comparison Table[l]shows the face-probabilities predicted by the Simp- 
son model (explained in Section [2]). This model fits neither of the two experi- 
ments. In comparison to both experiments, it clearly overpredicts the proba- 
bilities of the smallest faces and underestimates the probabilities of the largest 
faces. A much more accurate description is offered by the Gibbs models in Ta- 
ble[l] whose free parameter has been fitted to experiment I and II, respectively. 
This new model is now explained in Section |4j 
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4 Gibbs distribution 

Following independent ideas of Riemer (1991) and Obreschkow (2006), this sec- 
tion uses Gibbs distributions to model the face-probabilities of tossed cuboids. 
Gibbs distributions are probability distributions that are commonly used in 
many fields of probability theory, mathematical statistics, as well as statistical 
mechanics, from where they originate. The philosophy of this paper is to adopt 
Gibbs distributions in a heuristic way, that is without deriving them from a set 
physical assumptions. The model then gains its validity a posteriori though 
verification against experimental data - a common approach in statistics. 

A Gibbs distribution can be summarized as follows: consider a system with 
k states, where each state i = 1, ...,k has a positive energy E+. If the Gibbs 
theory applies, the system is found in state i with probability 

Pi (p) = Z(l3)- 1 e^>(-pE i ), (1) 

where j3 is a positive parameter, called inverse temperature (because it is 
proportional to T _1 in thermodynamics), and Z(/3) = J^- exp(— (3Ei) is a 
normalization factor, called the partition function. The parameter /3 controls 
the character of the Gibbs distribution: if /3 = the distribution is uniform 
with equal probabilities for all states i £ {1, k}; as /3 — > oo the distribution 
becomes peaked with the minimal energy state (s) having a probability equal 
to 1; for any intermediate /3 g (0,oo), the probability of a state increases 
monotonically with decreasing energy. 

In modeling the tossing experiments, the states are the faces that end up 
lying on top, i.e. the cuboid is said to be in state i if it comes to rest with face i 
on top (thus k = 6) . The energy of state i is taken proportional to the potential 
energy, i.e. to the height hi of the center of gravity in state i. Note that in this 
way inhomogeneities in the mass distribution of the cuboid are accounted for, 
as illustrated in Section [6j If the cuboid is homogeneous, hi = si/2 where Si 
is the vertical side-length of the cuboid in state i. To eliminate physical units, 
the energy Ei oc hi is normalized to the half-diagonal, 

Ei = ^ = ^ . (2) 

(Eli^) 1/3 (E-=i*i) 1/3 

Given this definition of the energies E il j3 is the only free parameter in the 
Gibbs distribution of Eq. ([I]). This parameter can be fitted to experimental 
data, for example using a maximum likelihood estimation (MLE). This method 
consists in maximizing the likelihood function 

6 6 

L(f3) ee H\pi(P)] ni = Z(P)~ N l[exp(-PEini), (3) 

i=i i=i 

where rii is the number of observations of state i and N = JZj—i n i- I n prac- 
tice, L(/3) can easily be maximized by minimizing In Z(f3) + /3J^ =1 Eifi. Dif- 
ferent values of /3 can maximize L(j3) in different experimental conditions. As 
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explained above, small values of j3 are expected in experimental conditions 
where all six faces appear frequently, while higher values of ft are expected if 
the smallest faces appear rarely. Explicitly, in the control experiment I (tossing 
with a dice cup) smaller values of (3 are expected, than in the control exper- 
iment II (free fall from 1 m height). In fact, the MLE yields (3 — 4.90 and 
j3 = 10.2 for experiments I and II, respectively. The corresponding probabil- 
ities of the Gibbs model are displayed in the bottom rows of Table [T] These 
probabilities (e.g. 43.0% for faces 3 and 4 in experiment I) lie often between 
the measured frequencies of the corresponding faces (e.g. 43.5% and 42.5%), 
thus suggesting that the model sufficiently describes the data. An explicit \ 2 
goodness-of-fit test shows that the predictions of the Gibbs model are indeed 
statistically consistent with the experimental data. This test will be used again 
and explained in more detail in the following section. 



5 Two classical experiments revisited 

Section [4] revealed that the Gibbs model offers a good approximation of the 
data gathered in the control experiment. The control experiment was based 
on a single cuboid with three different side-lengths. This section confronts the 
Gibbs model with other experimental data drawn from the literature. The 
main purpose of this comparison is to test whether the Gibbs model with a 
constant parameter f3 can describe a variety of differently shaped cuboids, 
tossed in similar experimental conditions. 

The experiments considered here are summarized in Table [2j They were 
performed by Budden (1980) and Heilbronner (1985), respectively. Both au- 
thors used families of ray-cuboids, i.e. cuboids with equal side-lengths s x in 
two orthogonal directions and a different side-length s y in the third direc- 
tion. Both authors experimented with a family of to cuboids j = l,...,m 
(to =15 for Budden and to = 7 for Heilbronner) with identical side- lengths 
Sx (s x = 15 mm for Budden and s x = 25 mm for Heilbronner), and vary- 
ing side-lengths s y j. Budden used cuboids "cut from a mild steel bar whose 
cross-section was a square of side 15 mm. These were distributed to a class 
of boys who tossed and rolled them while recording the results." By contrast, 
Heilbronner used cuboids from polyvinylchloride of density 1,500 kg m~ 3 . 
They were tossed "in the usual manner, i.e. rolled manually or from a shaker 
on cloth covered surfaces as well as on linoleum in a ratio of approximately one 
to one for each set of dice." The vast differences in the material and tossing 
techniques between Budden and Heilbronner suggest that these two datasets 
are described by different values (3 in the Gibbs model. However, the question 
to be investigated is whether within each dataset (Budden or Heilbronner) the 
to cuboids can be described by a constant parameter /3. 

To model the face-probabilities of a single ray-cuboid, the formalism of 
Section [4] can be simplified, since ray-cuboids only exhibit two macro-states: 
the ra-state, showing one of the two square faces, and the a;y-state showing 
one of the four rectangular faces. Given N tosses and n xx observations of the 
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xa;-state, the corresponding frequency is f xx = n xx /N . To model the face- 
probability p xx , the Gibbs model of Eq. ([I]) can be rewritten as 

Pxx = Z{P)- 1 eM-PE y ) (4) 

where E x = s x /(s x s x s y ) 1 ^ 3 , E y = s y /{s x s x s y ) x / 3 , and Z{0) = 2exp(-/3E x ) + 
exp(—(3E y ). Using Eq. ^ and p xy = 1 — p xx the likelihood function for (3 
becomes 

L(J3) =p xx (P) n **(l-p xx (P)) N - n **. (5) 

Given a set of m differently sized rrzy-cuboids j — 1, ...,m (with respec- 
tive variables p xx ,j, Nj, n XXt j, etc.), the best fitting constant parameter /? for 
all cuboids in the same dataset can be obtained by maximizing the global 
likelihood function 

m m 

m = n l m = iij , «jW b " j ( i -^w)^- ( 6 ) 

In practice, £(/3) can be maximized more easily by minimizing the function 
EjLitWj ln ^(/ ? ) + P(n xx ,jEy d + (1 - r^xj)^)]. This results in /3 = 4.46 
for the experiments performed by Budden and (3 — 3.53 for those performed 
by Heilbronner. The corresponding face-probabilities predicted by the Gibbs 
model are listed in the last column of Tableland plotted against s y /s x in 





Sa: [mm] 


s y [mm] 


N 




fm. [%] 


Pxx [%] 




15 


7.1 


332 


304 


91.6 


91.0 




15 


9.5 


840 


620 


73.8 


77.0 




15 


11.2 


799 


438 


54.8 


63.5 




15 


12.15 


740 


367 


49.6 


55.4 




15 


13.95 


516 


206 


39.9 


40.8 




15 


14.5 


530 


204 


38.5 


36.8 


o 


15 


17.4 


1011 


150 


14.8 


20.2 


-o 
-o 
~ 


15 


18.45 


532 


82 


15.4 


16.1 


pq 


15 


21.6 


654 


34 


5.2 


8.1 




15 


23.25 


606 


24 


4.0 


5.7 




15 


24 


702 


12 


1.7 


4.8 




15 


25.6 


609 


19 


3.1 


3.5 




15 


28 


680 


6 


0.9 


2.1 




15 


31.75 


275 


2 


0.7 


1.0 




15 


39.7 


503 


3 


0.6 


0.2 




25 


5 


2145 


2089 


97.4 


98.4 


CD 


25 


10 


2184 


1929 


88.3 


89.8 


d 
d 


25 


15 


2103 


1559 


74.1 


72.7 


o 
u 


25 


20 


2238 


1244 


55.6 


51.7 




25 


30 


2202 


421 


19.1 


20.5 




25 


35 


2259 


239 


10.6 


12.4 




25 


40 


2250 


162 


7.2 


7.6 



Table 2 Experimental datasets obtained by Budden and Heilbronner using different xxy- 
cuboids. The last column is the prediction of the Gibbs model using = AA6 (Budden) and 
P = 3.53 (Heilbronner). 
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Fig. 1 Measured frequencies f xx and fitted Gibbs probabilities p xx with a constant /3, as 
a function of the side-ratio s y /s x . The experimental values are those listed in Table [2] 
Vertical error bars represent standard deviations of f xx approximated as \J f xx /n xx , and 
horizontal error bars are standard deviations associated with 5% manufacturing errors for 
the side-lengths. 



Figure [TJ Qualitatively, there seems to be a good agreement between the ex- 
perimental data and the model. 

The rest of this section investigates whether the data is indeed statistically 
consistent - in a quantitative way - with the Gibbs model using a constant /3 
per dataset. To do so, a \ 2 test inspired by Gibbons and Chakraborti (2003) 
is used. The expected number of appearances of the rrcc-state is Njp xx j with 
a variance of Njp xx j. Therefore, the variance between experimental counts 
and model-prediction, normalized to the model variance, reads (Njp xx j — 
n xx.j) 2 /(Njp XX j). Applying an analogous reasoning to the xy-state and sum- 
ming over all the different cuboids, yields 

(Njp xxJ ~ n xx ^) 2 + (Njjl - p xxJ ) - (Nj - n XXJ )) 2 ~ 
N jPxx, 3 Nj(l-p xx j) \' 

If x 2 /m < 1, then the Gibbs model with a constant (3 per dataset fully de- 
scribes the experimental data; if x 2 / m > 1j the experimental data is not 
sufficiently matched by the model. Explicit calculations yield x 2 / m = 6.2 
(Budden) and \ 2 jm ~ 6.6 (Heilbronner), thus rejecting the hypothesis of a 
Gibbs model with a constant f3 per dataset at 2.5 standard deviations. In 
other words, this hypothesis seems to be rejected with a certainty of nearly 
99%. However, this \ 2 test ignores experimental uncertainties of various kinds, 
which shall now be addressed approximately. 

Potentially, there are various sources of systematic uncertainties in the ex- 
perimental data. For example, the tossing techniques might have been different 
for every xccy-cuboid. This is particularly plausible in Budden's experiment, 
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where different cuboids were tossed by different children. Further, the cuboids 
are not perfect due to material and manufacturing errors. To estimate the effect 
of manufacturing errors on the consistency between the data and the Gibbs 
model, an extended \ 2 test is performed, which explicitly accounts for uncer- 
tainties in the side-lengths s x and s y _j. It is assumed that these side-lengths 
are only known up to Gaussian errors with standard deviations es x and es y ,j. 
Hence, e represents the relative uncertainty of the side-lengths. The hypoth- 
esis that the m empirical f xx ,j belong to Gibbs distributions with a constant 
P is then tested using a parametric bootstrap test based on 999 independent 
iterations. In each iteration the following steps are executed. 

— For every cuboid j = I, m, 

— chose side-lengths s* = G(s x ,es x ) and s* ■ = G(s y j,es y j), where 
G{x, a) denotes a random number from a normal distribution with 
mean x and standard deviation er, 

— calculate the probability p* xx j of the xx-st&te using eq. Q with the 
original (3 (4.46 for Budden, 3.53 for Heilbronner) and the new side- 
lengths s* and s* p 

— simulate Nj tossing events, in which the irz-state appears with proba- 
bility p* xx j , and count the number n* xx j , 

— calculate the corresponding frequencies f xx j — n XX j/Nj. 

— Use the m values of f* x - to estimate the best parameter j3* via eq. (|3j). 

— Calculate the new probabilities p X x,j using eq. Q with j3* , s x and s y j. 

— Calculate the value \ 2 using eq. (J7]) with the probabilities p xx j. 

If the original % 2 is large in comparison to the 999 simulated values of \ 2 > 
the hypothesis of a constant /3 must be rejected. 

However, the p- values listed in Table[3]show that already values around e = 
0.05, i.e. manufacturing errors of 5%, make the experimental data compatible 
with the hypothesis of a constant /3 for all m cuboids. Measurements of the 
masses of machine-manufactured wood cuboids similar to those of Budden 
revealed mass deviations around 7% between 'identical' cuboids, roughly in 
line with side-length variations of 5%. 

In summary, the Gibbs model with a constant ft for all the cuboids in a 
dataset (Budden or Heilbronner) is consistent with the experimental data as 
long as plausible manufacturing errors are accounted for. Figuratively speak- 
ing, the data points in Figure [T] are consistent with the models, as long as the 
horizontal error bars are included. 



e 


Budden 


Heilbronner 


0.03 


0.000 


0.003 


0.04 


0.006 


0.021 


0.05 


0.067 


0.090 


0.06 


0.187 


0.206 



Table 3 p- values of x 2 in the simulated \ 2 -distribution to test the hypothesis of the Gibbs 
model with a constant /3 per dataset (Budden or Heilbronner). 



Cuboidal Dice and Gibbs Distributions 



9 



6 Extension to non-cuboidal dice 

As shown so far, the Gibbs model fully describes the face-probabilities of tossed 
cuboids within the uncertainties of currently available experimental data. This 
motivates the idea that the Gibbs model could be extended to more complex 
dice geometries and inhomogeneous cuboids. A full investigation of this idea 
lies beyond the scope of this paper, but to provide an illustration the U-shaped 
die shown in Figure [2] is considered. Two experimental runs were performed 
with this die. In experiment I, the die was tossed N = 1, 950 times onto a hard 
surface; in experiment II it was dropped N = 150 times onto a wool carpet. 
The measured frequencies of the different faces are listed in Table |4j Unlike the 
cuboid, the U-shaped die has no symmetry between the faces 3 and 4. However, 
the Gibbs model as given in Eq. can still be applied using the heights hi of 
the center of gravity listed in Table [3] To calculate the corresponding energies 
Ei the hi are normalized to the half-diagonal of 16.45 mm. The maximum 
likelihood method yields f3 = 5.11 (experiment I) and j3 = 8.41 (experiment 
II), respectively. The higher j3 of the second experiment is clearly related to 
the wool carpet's softness, which tends to destabilize positions with a high 
center of gravity, thus making the probabilities more skewed towards the most 
stable positions. 

The probabilities of the Gibbs model are consistent with the data in terms 
of the x 2 test discussed in Section [5j hence demonstrating that the Gibbs 
model extends to non-cuboidal dice. 




Fig. 2 Image of the U-shaped die. The digits are the indices of the visible faces. Face 1 
(opposite face 6) is hidden on the left, face 2 (opposite face 5) is hidden at the back, and 
face 3 (opposite face 4) is hidden at the bottom. 
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Face i 


1 


2 


3 


4 


5 


6 


Heights of center of gravity hi [mm] 


10.0 


11.5 


7.61 


5.39 


11.5 


10.0 


fi experiment I (N = 1, 950) [%] 
fi experiment II (N = 150) [%] 


10.6 
4.7 


6.9 
2.0 


23.9 
28.0 


42.5 
57.3 


6.8 
1.3 


9.3 
6.7 


Pi Gibbs model (/3 = 5.11) [%] 
Pi Gibbs model (/3 = 8.41) [%] 


10.4 
5.9 


7.3 
3.3 


21.9 
20.0 


43.6 
62.2 


6.5 
2.7 


10.4 
5.9 



Table 4 Results of tossing the U-shaped die shown in Figure |2 Note that experiment II 
has a very small number TV, thus very large statistical uncertainties on the values 



7 Summary 

This paper uncovered that the face-probabilities of a tossed cuboid are well de- 
scribed by the Gibbs model defined via Eqs. and p]). These face-probabilities 
depend heavily on the tossing conditions - an effect that can be accounted for 
by the Gibbs model by adjusting the free parameter f3. Good fits of /3 can 
be obtained via the maximum likelihood method of Eq. ([3|. Typical values 
of f3 range between 3 and 10. If differently shaped cuboids are all tossed us- 
ing similar conditions (material, technique, etc.), then the face-probabilities 
of all these cuboids can be well approximated using a constant parameter /3, 
estimated via the global maximum likelihood method of Eq. 

Acknowledgements The authors thank Robert Allin for valuable discussions about an 
earlier version of this paper. D.O. acknowledges the discussions with Nick Jones. 



References 

1. Budden, F., "Throwing non-cubical dice", The Mathematical Gazette, 64(429), 196-198 
(1980) 

2. Gibbons, J. D. and Chakraborti, S. "Nonparametric Statistical Inference", 4th Edition, 
Marcel Dekker, New York, Basel (2003) 

3. Heilbronner, E. "Crooked dice", J. Recreational Math. 17, 177 (1985) 

4. Hyskova, M. and Kalousova, A. and Saxl, I. llEarly history of geometric probability and 
stereology", Image Anal Stereol 31, 1-16 (2012) 

5. Newton, I. and Whiteside D. T. (ed.) "The Mathematical Papers of Isaac Newton", vol. 
I. 1664-1666, pages 60-61 Cambridge University Press, Cambridge (1967) 

6. Obreschkow, D. "Broken symmetry and the magic of irregular dice", Online project: 
http://www. quantumholism. com/cuboid (2006) 

7. Riemer, W. "Stochastische Probleme aus elementarer Sicht", Bibliographisches Institut, 
Mannheim, Wien, Ziirichm (1991) 

8. Simpson, T. "The Nature and Laws of Chance" Cave, London (1740) 

9. Singmaster D. "Theoretical probabilities for a cubical die" The Math. Gazette. 65, 
208-210 (1981) 



